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Abstract 


A  new  technique  for  approximating  the  scattering  of  laser  light  by  biological  cells  is  reported. 
This  technique  is  based  on  a  three-dimensional  scalar  wave  equation  approximation  of  the  full 
Maxwell’s  field  equations  for  the  electromagnetic  field.  This  scalar  wave  equation  for 
describing  the  light  scattering  patterns  of  cells  containing  arbitrary  morphological  structure 
(e.g.,  various  organelles)  is  solved  numerically  using  a  spectral  method.  The  accuracy  of  the 
spectral  numerical  method  is  verified  by  comparison  with  solutions  obtained  from  linear 
perturbation  theory  and  Mie  theory.  Comparison  with  Mie  theory  shows  that  the 
three-dimensional  scalar  wave  equation  is  a  good  approximation  to  the  full  Maxwell’s  field 
equations  for  light  scattering  up  to  moderate  forward  scattering  angles  (i.e.,  for  scattering 
angles  less  than  about  35°).  The  approximate  technique  used  here  is  capable  of  correctly 
predicting  the  scattered  intensity  patterns  from  biological  cells  over  a  dynamic  range  spanning 
six  orders  of  magnitude.  The  new  technique  can  be  applied  to  calculate  the  light  scattering 
either  from  an  individual  biological  cell  and  from  a  sample  containing  an  ensemble  of  such 
biological  cells.  The  scattering  intensity  patterns  predicted  using  the  new  technique  can 
potentially  be  applied  to  diagnose  the  size  and  internal  structure  of  biological  cells,  making  it  a 
valuable  interpretative  tool  in  flow  cytometry  (e.g.,  in  the  detection  of  rare  event  cells  such  as 
those  resulting  from  a  biological  warfare  agent  attack,  or  for  the  rapid  noninvasive  optical 
assessment  of  tissue  pathology  in  the  detection  of  cancerous  cells). 
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Resume 


On  documente  ici  une  nouvelle  technique  faisant  une  approximation  de  la  diffusion  de  la 
lumiere  laser  au  moyen  de  cellules  biologiques.  Cette  technique  est  basee  sur 
P approximation  d’une  equation  d’onde  scalaire  tri  dimensionnelle  des  equations  completes 
de  Maxwell  pour  le  champ  electromagnetique.  Cette  equation  d’onde  scalaire  visant  a 
decrire  les  modeles  de  diffusion  de  la  lumiere  de  cellules  contenant  une  structure 
morphologique  (par  ex. :  organelles  variees)  est  resolue  numeriquement  au  moyen  d’une 
methode  spectrale.  L’exactitude  de  cette  methode  numerique  spectrale  est  verifiee  en 
comparant  avec  des  solutions  obtenues  a  partir  de  la  theorie  de  la  perturbation  lineaire  et  de 
la  theorie  Mie.  La  comparaison  avec  la  theorie  Mie  indique  que  l’equation  d’onde  scalaire 
tri  dimensionnelle  est  une  bonne  approximation  des  equations  completes  de  Maxwell,  en  ce 
qui  conceme  la  diffusion  de  la  lumiere,  tant  que  ces  angles  diffusent  moderement  vers 
l’avant  (angles  de  diffusion  de  moins  de  35°  a  peu  pres).  La  technique  appropriee  utilisee  ici 
est  capable  de  predire  correctement  les  modeles  d’intensite  diffuses  a  partir  des  cellules 
biologiques  a  une  echelle  dynamique  comprenant  six  ordres  de  magnitude.  La  nouvelle 
technique  peut  etre  appliquee  ici  pour  calculer  la  diffusion  de  la  lumiere  a  partir  aussi  bien 
d’une  cellule  biologique  individuelle  que  d’un  echantillon  contenant  un  ensemble  de  ces 
cellules  biologiques.  Les  modeles  d’intensite  de  diffusion,  predits  en  utilisant  la  nouvelle 
technique,  peuvent  etre  appliques  au  diagnostic  de  la  taille  et  de  la  structure  interne  des 
cellules  biologiques  et  devenir  un  outil  d’interpretation  precieux  en  cytometrie  de  flux  (par 
ex.  :  dans  la  detection  de  cellules  rares  telles  que  celles  resultant  d’une  attaque  d’agent  de 
guerre  biologique  ou  bien  pour  1’ evaluation  optique  rapide  non  effractive  de  la  pathologie 
des  tissus,  en  ce  qui  conceme  la  detection  de  cellules  cancereuses. 
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Executive  summary 

Background 


A  high  priority  of  any  biological  warfare  (BW)  agent  defence  system  is  the  triggering/cueing 
of  the  system  in  the  presence  of  incoming  aerosol  clouds  that  may  contain  the  presence  of  BW 
agents,  followed  by  the  initiation  of  particulate  air  sampling  required  for  the  subsequent 
detection  and  identification  of  the  putative  agent  in  order  to  warn  the  user  of  the  presence  of  a 
BW  agent  cloud.  In  this  regard,  the  problem  of  the  detection  of  BW  agents  (e.g.,  bacteria, 
viruses)  based  on  some  small  set  of  quantifiable  properties  that  can  be  shown  to  differ  from 
their  immediate  background  clutter  is  of  great  interest.  Indeed,  this  problem  was  emphasized 
by  Brigadier  General  Doesberg  (Commander,  Joint  Program  Office  for  Biological  Detection, 
Alexandria,  VA)  who  in  his  keynote  address  to  a  1 996  Edgewood  Research,  Development,  and 
Engineering  Center  (ERDEC)  scientific  conference  emphasized  that  “the  major  deficiency  of 
our  biological  detection  system  is  the  lack  of  a  reliable  triggering/detection  device”. 

The  key  problem  that  arises  in  the  detection  of  a  BW  agent  is  to  distinguish  unambiguously  the 
putative  presence  of  a  pathogenic  bio-target  (BW  agent)  against  the  ambient  background 
clutter  (which  can  consist  of  inorganic  matter  such  as  sulfates  and  nitrates  and  live 
non-pathogenic  matter  such  as  pollen,  fungi  and  bacteria).  Most  BW  agent  detection  systems 
involve  drawing  an  ambient  air  sample  through  a  laser-illuminated  sampling  volume  where 
airborne  particles  scatter  light  and/or  undergo  particle  fluorescence.  Hence,  increasing  the 
present  understanding  of  the  interaction  of  laser  light  with  biological  cells  can  potentially 
improve  the  ability  of  current  BW  agent  detection  systems  to  discriminate  between  a  target 
BW  agent  and  environmental  interferents.  Understanding  the  relationship  between  measured 
scattered  light  properties  and  physiological  differences  at  the  cellular  level  is  fundamental  to 
providing  a  high  degree  of  selectivity  for  BW  agents  in  current  BW  detection  systems  and 
constitutes  the  main  challenge  for  the  modelling  of  light  scattering  from  biological  cells. 


Principal  Results 

In  this  report,  we  develop  a  simplified  mathematical  model  for  the  interaction  of  laser  light 
with  biological  cells.  In  this  approach,  the  full  set  of  Maxwell’s  field  equations  is  reduced  to  a 
wave  equation  for  the  electric  field.  This  wave  equation  is  solved  numerically  using  a  spectral 
method.  This  method  provides  a  computationally  efficient  numerical  scheme  that  can  be  used 
to  predict  the  light  scattering  intensity  pattern  from  either  a  single  biological  cell  or  a  sample 
containing  many  biological  cells. 

The  accuracy  of  the  spectral  numerical  method  is  verified  by  comparison  with  solutions 
obtained  from  linear  perturbation  theory  and  Mie  theory.  Comparison  with  Mie  theory  shows 
that  the  three-dimensional  scalar  wave  equation  is  a  good  approximation  to  the  full  Maxwell’s 
field  equations  for  light  scattering  up  to  moderate  forward  scattering  angles  (i.e.,  for  scattering 
angles  less  than  about  35s).  The  approximate  technique  used  here  is  capable  of  correctly 
predicting  the  scattered  intensity  patterns  from  biological  cells  over  a  dynamic  range  spanning 
six  orders  of  magnitude.  The  new  technique  can  be  applied  to  calculate  the  light  scattering 
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either  from  an  individual  biological  cell  and  from  a  sample  containing  an  ensemble  of  such 
biological  cells.  The  scattering  intensity  patterns  predicted  using  the  new  technique  can 
potentially  be  applied  to  diagnose  the  size  and  internal  structure  of  biological  cells,  making  it  a 
valuable  interpretative  tool  in  flow  cytometry  (e.g.,  in  the  detection  of  rare  event  cells  such  as 
those  resulting  from  a  biological  warfare  agent  attack,  or  for  the  rapid  noninvasive  optical 
assessment  of  tissue  pathology  in  the  detection  of  cancerous  cells). 


Significance  and  Future  Plans 

To  improve  the  performance  of  BW  agent  detection  systems,  it  is  necessary  to  understand  the 
relationship  between  optical  and  biological  properties  of  BW  agents  (bacteria,  viruses, 
ridkettsiae).  A  model  designed  to  provide  insight  into  the  relationship  between  the  scattering 
properties  of  a  biological  cell  and  its  morphology  has  been  developed  in  this  study.  The  model 
can  potentially  be  used  to  provide  deeper  insights  into  the  characteristics  of  laser  light 
scattering  from  biological  cells.  In  consequence,  this  fundamental  understanding  may  lead  to 
an  improved  ability  to  discriminate  biological  agents  from  other  harmless  biological  and 
non-biological  materials  present  in  the  ambient  environment. 

The  model  developed  here  provides  a  capability  for  interpreting  laser  light  scattering  patterns 
from  BW  agents,  and  in  particular  for  assessing  cellular  morphology  from  scattering 
measurements.  This  will  lead  to  the  development  of  novel  analysis  methods  for  flow  cytometry 
which  in  turn  can  lead  to  improvements  in  BW  agent  detection.  Future  work  should  focus  on 
incorporating  a  detailed  numerical  model  for  prediction  of  light  scattering  from  biological 
cells  with  a  new  generation  of  micro-cytometers  (viz.,  flow  cytometry  undertaken  on  a 
microchip)  that  would  allow  a  more  precise  interpretation  of  the  data  obtained  from  such  a 
device  (and,  particularly,  in  relation  to  the  extraction  of  morphological  parameters  from  BW 
agent  cells  based  on  scattering  measurements). 


Shao,  Y.  and  Yee,  E.  (2004).  An  Approximate  Method  for  Modelling  Laser  Light  Scattering 
from  Biological  Cells.  (DRDC  Suffield  TR  2004-187).  Defence  R&D  Canada  -  Suffield. 
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Sommaire 


Contexte 

La  grande  priorite  d’un  systeme  de  defense  contre  un  agent  de  guerre  biologique  (BW)  est 
d’etre  en  mesure  de  se  declencher  ou  de  signaler  l’arrivee  de  nuages  aerosol  pouvant  contenir 
des  agents  BW,  puis  d’amorcer  l’echantillonnage  d’air  particulier  requis  pour  detecter  et 
identifier  ulterieurement  l’agent  putatif,  ceci  visant  a  prevenir  l’utilisateur  de  la  presence  d’un 
nuage  d’agents  BW.  A  cet  egard,  il  devient  tres  interessant  d’etudier  le  probleme  de  la 
detection  des  agents  BW  (par  ex.  :  bacteries,  virus)  en  se  basant  sur  un  certain  petit  ensemble 
de  proprietes  quantifiables  pouvant  demontrer  qu’un  agent  differe  du  fouillis  de  son 
environnement  immediat.  Ce  probleme  a  ete  certainement  souligne  par  le  brigadier  general 
Doesberg  (Commandant  du  Bureau  commun  du  programme  de  detection  biologique, 
Alexandria,  VA).  Dans  son  discours  liminaire  a  la  conference  scientifique  de  Edgewood 
Research,  Development,  and  Engineering  Center  (ERDEC),  en  1996,  il  a  souligne  que 
«  1’ inconvenient  majeur  de  notre  systeme  de  detection  biologique  est  le  manque  d’appareil  de 
declenchement  et  de  detection  fiable  ».  Le  probleme  cle  de  la  detection  d’un  agent  BW  est  la 
difficult^  de  distinguer  sans  ambigui'te  la  presence  putative  d’une  cible  biologique 
pathogenique  (agent  BW)  dans  le  fouillis  ambiant  de  son  environnement  (pouvant  consister  en 
matiere  non  organique  telle  que  les  sulfates  et  les  nitrates  et  de  matiere  vivante  non 
pathogenique  telle  que  le  pollen,  champignons  microscopiques  et  bacteries).  La  plupart  des 
systemes  de  detection  d’agents  BW  consistent  a  retirer  un  echantillon  d’air  ambiant  a  l’aide 
d’un  volume  d’echantillon  illumine  au  laser  dans  lequel  les  particules  aerosol  diffusent  la 
lumiere  et  /  ou  la  particule  devient  fluorescente.  Ainsi,  une  meilleure  comprehension  de 
1’ interaction  entre  la  lumiere  laser  et  les  cellules  biologiques  peut  sans  doute  ameliorer  la 
capacite  des  systemes  actuels  de  detection  d’agents  BW  a  discriminer  entre  un  agent  BW  cible 
et  ses  interferents  environnementaux.  Il  est  fondamental  de  comprendre  la  relation  entre  les 
proprietes  de  la  lumiere  diffusee  mesuree  et  les  differences  physiologiques  au  niveau 
cellulaire,  pour  foumir  un  haut  niveau  de  selectivity,  en  ce  qui  conceme  les  agents  BW  dans 
les  systemes  actuels  de  detection  BW  et  constitue  le  defi  principal  dans  le  domaine  de  la 
modelisation  de  la  diffusion  de  la  lumiere.  a  partir  de  cellules  biologiques. 

Resultats  principaux 

Dans  ce  rapport,  nous  developpons  un  modele  mathematique  simplifie  de  l’interaction  de  la 
lumiere  laser  avec  les  cellules  biologiques.  Cette  methode  reduit  l’ensemble  complet  des 
equations  de  Maxwell  a  une  equation  d’onde  pour  le  champ  electrique.  L’equation  d’onde  est 
resolue  numeriquement  en  utilisant  une  methode  spectrale.  Cette  methode  foumit  un  schema 
de  calcul  numerique  efficace  pouvant  etre  utilise  pour  predire  le  modele  de  l’intensite  de  la 
diffusion  de  la  lumiere  a  partir  soit  d’une  seule  cellule  biologique  soit  d’un  echantillon 
contenant  beaucoup  de  cellules  biologiques. 

L’exactitude  de  cette  methode  numerique  spectrale  est  verifiee  en  comparant  avec  des 
solutions  obtenues  a  partir  de  la  theorie  de  la  perturbation  lineaire  et  de  la  theorie  Mie.  La 
comparaison  avec  la  theorie  Mie  indique  que  l’equation  d’onde  scalaire  tri  dimensionnelle  est 
une  bonne  approximation  des  equations  completes  de  Maxwell,  en  ce  qui  conceme  la 
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diffusion  de  la  lumiere,  tant  que  ces  angles  diffusent  moderement  vers  l’avant  (angles  de 
diffusion  de  moins  de  35°  a  peu  pres).  La  technique  appropriee  utilisee  ici  est  capable  de 
predire  correctement  les  modeles  d’intensite  diffusee  a  partir  des  cellules  biologiques  a  une 
echelle  dynamique  comprenant  six  ordres  de  magnitude.  La  nouvelle  technique  peut  etre 
appliquee  ici  pour  calculer  la  diffusion  de  la  lumiere  a  partir  aussi  bien  d  une  cellule 
biologique  individuelle  que  d’un  echantillon  contenant  un  ensemble  de  ces  cellules 
biologiques.  Les  modeles  d’intensite  de  diffusion,  predits  en  utilisant  la  nouvelle  technique, 
peuvent  etre  appliques  au  diagnostic  de  la  taille  et  de  la  structure  interne  des  cellules 
biologiques,  ce  qui  les  rend  un  outil  d’interpretation  precieux  en  cytometrie  de  flux  (par  ex.  : 
dans  la  detection  de  cellules  rares  telles  que  celles  resultant  d’une  attaque  d’ agent  de  guerre 
biologique  ou  bien  encore  pour  revaluation  optique  rapide  non  effractive  de  la  pathologie  des 
tissus,  en  ce  qui  conceme  la  detection  de  cellules  cancereuses. 

La  portee  des  resultats  et  les  plans  futurs 

II  est  necessaire  de  comprendre  la  relation  entre  les  proprietes  optiques  et  biologiques  des 
agents  BW  (bacterie,  virus,  rickettsie)  pour  ameliorer  la  performance  des  systemes  de 
detection  d’agents  BW.  Cette  etude  a  developpe  un  modele  con?u  pour  procurer  une 
connaissance  plus  approfondie  de  la  relation  entre  les  proprietes  de  diffusion  d  une  cellule 
biologique  et  sa  morphologie.  II  est  possible  que  ce  modele  soit  utilise  pour  acquerir  une 
connaissance  plus  profonde  des  caracteristiques  de  la  diffusion  de  la  lumiere  laser  a  partir  de 
cellules  biologiques.  Cette  connaissance  fondamentale  peut,  par  consequent,  aboutir  a 
l’amelioration  de  la  capacite  de  discriminer  les  agents  biologiques  de  ceux  qui  sont  moins 
dangereux  et  des  materiaux  non  biologiques  qui  sont  presents  dans  l’environnement  ambiant. 

Le  modele  developpe  ici  donne  la  capacite  d’interpreter  les  modeles  de  diffusion  de  lumiere 
laser  a  partir  d’agents  BW.  II  permet  surtout  d’evaluer  la  moiphologie  cellulaire  a  partir  des 
mesures  prises  pendant  la  diffusion.  Ceci  amenera  a  la  mise  au  point  de  methodes  d  analyse 
nouvelles  pour  la  cytometrie  de  flux  qui  a  son  tour  peut  aboutir  a  l’amelioration  de  la 
detection  d’agent  BW.  Les  travaux  futurs  devraient  se  concentrer  sur  1’ incorporation  de 
modeles  numeriques  detailles  pour  la  prediction  de  la  diffusion  de  la  lumiere  a  partir  de 
cellules  biologiques  avec  une  nouvelle  generation  de  micro-cytometres  (viz.,  cytometrie  de 
flux  entreprise  sur  un  microcircuit  integre)  qui  permettrait  une  interpretation  plus  precise  des 
donnees  obtenues  a  partir  d’un  tel  appareil  (et  surtout  celles  en  rapport  avec  1  extraction  de 
parametres  morphologiques  a  partir  de  cellules  d’agents  BW,  basee  sur  les  mesures  prises 
durant  la  diffusion). 
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Introduction 


Lasers  have  several  unique  properties,  such  as  monochromaticity,  high  power,  short  pulse 
duration,  and  coherence,  which  are  particularly  useful  for  biological  agent  detection  and 
bio-medical  applications.  Laser  medicine  is  now  a  rapidly  growing  research  field.  The 
increasing  use  of  laser  light  for  detection  of  biological  warfare  (BW)  agents  and  for  diagnostic 
and  therapeutic  medicine  has  created  a  need  to  understand  how  laser  light  interacts  with 
biological  cells  and  propagates  through  biological  tissue.  This  understanding  would  result  in 
improved  BW  agent  detection  by  providing  a  quantitative  interpretation  of  the  signals  that  are 
measured,  and  in  bio-medical  applications  would  enable  the  quantitative  analysis  of  diagnostic 
measurements  leading  potentially  to  an  optimal  development  of  therapeutic  techniques. 

There  are  five  main  categories  [1,  2,  3]  of  photo-physical  processes  in  the  interaction  of  laser 
light  with  biological  tissues.  These  are  non-destructive  photo,  photo-chemical,  photo-thermal, 
photo-ablative  and  photo-mechanical  interactions.  Non-destructive  photo  interactions  mainly 
concern  laser  light  propagation  and  scattering  in  biological  tissue,  which  can  be  used  both  for 
imaging  and  diagnostics.  Photo-chemical  interactions  involve  the  absorption  of  light  by 
specific  molecules  that  are  either  present  in,  or  added  to,  a  biological  tissue  sample.  Such 
interactions  are  the  basis  for  photo-dynamic  therapy.  Photo-thermal  interactions  are  those 
where  the  observed  biological  effect  is  due  to  the  deposition  of  heat  in  the  tissue.  Most  current 
laser  surgery,  such  as  welding  and  coagulation,  for  example,  falls  into  this  category. 
Photo-ablative  interactions  can  occur  when  photons  have  sufficient  energy  to  cause  the 
dissociation  of  biopolymers  with  a  subsequent  desorption  of  the  resulting  fragments.  The 
threshold  power  density  for  this  effect  to  occur  is  approximately  108  W  cm-2  for  10  ns  laser 
pulses.  Photo-mechanical  interactions  occur  at  fluence  rates  of  approximately  1010  W  cm-2 
for  nanosecond  laser  pulses  and  1012  W  cm-2  for  picosecond  pulses.  When  these  pulses 
illuminate  the  biological  tissue,  the  dielectric  of  the  tissue  experiences  a  breakdown,  and  a 
small  volume  plasma  is  produced.  The  expansion  of  this  plasma  creates  a  shock  wave  which 
can  mechanically  rupture  the  tissue.  These  last  two  types  of  interaction  are  complex  in  terms 
of  their  occurrence  thresholds  and  concomitant  nonlinear  effects. 

In  this  report,  we  will  focus  on  non-destructive  photo  interactions  of  laser  light  with  biological 
cells,  where  the  optical  properties  of  the  biological  cell  are  time  invariant  and  are  independent 
of  the  probing  light  field.  Optical  properties  of  biological  tissue  and  cells  are  very  important 
for  a  wide  range  of  studies  ranging  from  detection/imaging  and  diagnostic  applications,  such 
as  cytometry  [4],  confocal  [5]  and  optical  coherence  tomography  [6]  imaging,  and 
fundamental  investigations  of  the  sensitivity  of  a  biological  cell  to  light  sources  [7].  When 
laser  light  enters  a  biological  sample,  it  can  be  scattered  and  absorbed.  The  relative  probability 
of  these  processes  in  a  given  biological  cell  or  tissue  sample  depends  on  the  laser  wavelength 
[8].  In  these  applications,  laser  light  from  the  blue  to  near-infrared  parts  of  the  spectrum  is 
predominantly  scattered  from  local  variations  in  the  refractive  index  between  different  parts  of 
the  biological  cell. 

Cytometry  is  the  measurement  of  physical  and/or  chemical  characteristics  of  biological  cells 
(e.g.,  bacteria,  viruses,  etc.).  Flow  cytometry  [9]  is  a  process  in  which  such  measurements  are 
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made  while  the  cells  or  particles  are  passing  through  a  sampling  volume  in  a  fluid  stream. 
During  the  past  30  years,  sophisticated  improvements  have  made  flow  cytometry  a  powerful 
and  valuable  tool  for  the  quantitative  analysis  of  individual  cells  or  other  biological  particles. 

A  flow  cytometer  is  an  instrument,  which  can  be  used  to  obtain  quantitative  information  based 
on  light  scattering  or  fluorescence  emission  caused  by  individual  cells  (or  particles)  as  they 
flow  rapidly  in  a  fluid  stream  in  front  of  a  light  source  [4],  The  components  of  a  flow 
cytometer  usually  include  a  light  source  (normally  a  laser);  a  sample  chamber  with  a  flow  cell 
and  a  sheath  fluid  stream;  a  photodetector  or  photomultiplier  tube  (PMT)  that  collects  light  and 
converts  it  to  an  electronic  signal;  a  signal  processing  system  that  converts  the  signal  from  an 
analog  form  to  a  digital  format;  and  a  computer  to  direct  operations,  store  the  collected  signals, 
and  display  data  [10].  With  the  combination  of  state-of-the-art  advances  in  computer  and  laser 
technology,  more  sophisticated  flow  cytometers  have  been  developed  to  obtain  objective  and 
precise  measurements  of  multiple  characteristic  parameters  of  an  individual  cell  at  one  time. 
Such  parameters  include  cell  size,  cell  shape,  and  cell  refractive  index.  Recently,  flow 
cytometers  [e.g.,  Los  Alamos  National  Laboratory  flow  cytometer  and  the  Becton  Dickenson 
Flow  Cytometer  (FACSCaliber)]  have  been  used  in  bio-detection  to  determine  if  particulates 
collected  are  biological  or  inorganic  in  origin. 

Understanding  relationships  between  the  measured  scattered  light  properties  and  physiological 
characteristics  at  the  cellular  level  in  the  biological  cell  is  fundamental  to  the  usefulness  of 
optical  diagnostics  (e.g.,  use  of  flow  cytometry  for  bio-detection)  and  constitutes  the  main 
challenge  for  the  modelling  of  laser  light  scattering  from  cells.  The  complexity  of  the 
scattering  media,  coupled  with  comparable  spatial  scales  of  inhomogeneity  in  the  dielectric 
constant  and  the  laser  wavelength,  limits  an  analytical  approach  to  the  scattering  problem. 
Several  numerical  procedures  have  been  adoped  to  deal  with  the  laser  light-cell  scattering 
problem  from  their  original  physical  and  engineering  applications. 

The  most  complete  numerical  method  is  based  on  the  finite-difference  time-domain  (FDTD) 
algorithm,  which  was  originally  proposed  [1 1]  as  the  method  to  solve  the  full  set  of  Maxwell’s 
equations.  Since  this  seminal  work,  the  numerical  solution  of  Maxwell’s  equations  using  the 
FDTD  method  has  been  applied  to  perform  electromagnetic  simulations  in  many  areas  in 
electrical  engineering  and  related  fields  such  as  microwave  engineering,  remote  sensing  and 
subsurface  sensing,  antenna  analysis  and  design,  radar  cross-section  analysis  and  design, 
stealth  technology,  wireless  communication  and  propagation,  and  electromagnetic 
compatibility  and  electromagnetic  interference  (EMC/EMI).  It  is  only  very  recently  that  some 
preliminary  effort  has  been  expended  on  the  application  of  electromagnetic  simulations  to 
address  bioengineering  and  biotechnology  problems.  In  a  recent  implementation  of  the  FDTD 
algorithm  to  the  light-cellular  scattering  problem  [12],  [13],  the  biological  cell  has  been 
modelled  through  the  inhomogeneous  refractive  index  of  its  different  structures  (organelles). 
The  FDTD  approach  requires  large  computational  resources  both  in  terms  of  memory  and  time 
because  the  FDTD  algorithm  solves  for  the  full  vector  electric  and  magnetic  fields  and 
resolves  their  evolution  on  the  time  scale  of  a  laser  light  period. 

The  Monte  Carlo  (MC)  method  is  a  technique  first  proposed  to  simulate  physical  processes 
using  a  stochastic  model  [14]  and  it  has  been  used  to  solve  a  variety  of  physical  problems. 
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Recently,  MC  simulations  of  photon  propagation  provide  a  flexible,  yet  rigorous,  approach 
toward  photon  propagation  in  turbid  tissues  [15].  It  should  be  pointed  out  that  the  MC 
simulations  are  rigorous,  but  necessarily  statistical  in  nature  and,  therefore,  require  significant 
computation  time  to  achieve  the  required  precision  and  resolution.  An  alternative  to  the  FDTD 
and  MC  methods  is  the  discrete  dipole  approximation  [16],  [17],  in  which  a  particle  is  divided 
into  subvolumes  that  are  assumed  to  behave  as  dipoles. 


In  summarizing,  the  FDTD  and  MC  methods  are  the  two  mainstream  numerical  methods  for 
modelling  laser  light  propagation  and  scattering  from  biological  cells  and  tissues.  The  MC 
method  efficiently  describes  multiple  scattering  from  cells.  However,  it  is  not  well  suited  to 
addressing  scattering  from  a  single  cell.  On  the  contrary,  the  FDTD  method  accurately 
describes  scattering  from  a  single  cell,  but  its  application  to  multiple  cell  scattering  results  in 
computational  requirements  that  increase  dramatically  with  an  increase  in  the  size  of  the 
computational  domain.  In  view  of  this,  the  objectives  of  this  report  are:  (1)  to  develop  an 
approximate  technique  that  can  be  used  to  predict  the  scattering  intensity  patterns  from 
samples  consisting  of  single  or  multiple  biological  cell(s);  and,  (2)  to  develop,  implement,  and 
validate  a  computer  simulation  code  based  on  this  technique  for  the  modelling  and  prediction 
of  laser  light  scattering  from  biological  cells. 


The  proposed  new  approximate  method  for  modelling  laser  light  scattering  from  biological 
cells  is  based  on  a  simplified  mathematical  model  of  light  propagation.  The  resulting  equations 
of  this  method  are  solved  numerically  using  a  spectral  algorithm.  The  results  of  numerical  and 
analytical  studies  of  light  scattering  from  two-dimensional  (2-D)  and  three-dimensional  (3-D) 
cell  geometries  are  presented  in  this  report.  In  the  2-D  cellular  scattering  studies,  the  relative 
change  in  the  angular  distribution  of  scattering  intensity  due  to  the  internal  cell  structure, 
shape,  size,  and  refractive  index  can  be  addressed  more  easily  than  in  three  dimensions  because 
the  computational  burden  and  memory  requirements  are  smaller.  Examples  of  multiple 
scattering  processes  (i.e.,  dependence  of  the  scattered  light  characteristics  on  the  number  of 
scatterers  and  density  of  cells  in  the  sampling  volume)  are  provided  for  2-D  cell  geometries. 


In  our  approximate  simulation  approach,  the  full  set  of  Maxwell’s  equations  is  reduced  to  a 
wave  equation  for  the  electric  field.  Further  approximations  include  the  elimination  of  high 
frequency  (o>o)  oscillations  (carrier  frequency  of  the  electromagnetic  wave)  by  enveloping  the 
field  amplitude  variations  that  occur  on  an  (0q  1  time  scale,  which  reduces  the  order  of  the  time 
derivative  in  the  wave  equation  (a  detailed  explanation  of  this  process  is  provided  later  in 
Section  2).  As  compared,  for  example,  to  numerically  solving  the  full  set  of  Maxwell’s 
equations,  we  have  reduced  the  number  of  equations  to  solve,  eliminated  the  need  for  temporal 
resolution  on  an  (Oq  1  time  scale,  and  by  applying  the  spectral  method  [18]  for  the  numerical 
solution  of  the  resulting  equation  have  enabled  the  scattering  intensity  distribution  to  be 
accurately  calculated  using  six  grid  points  per  laser  wavelength.  By  reducing  the  numerical 
requirements  to  obtain  a  solution,  we  have  developed  an  effective  algorithm,  which  could  be 
applied  to  studies  of  scattering  from  a  single  biological  cell  as  well  as  from  a  sample 
consisting  of  ensemble  of  many  biological  cells. 
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The  Computational  Model  _ _ 

Basic  Equations  of  the  Approximate  Model 

The  computational  model  used  for  the  scattering  problem  considered  in  this  report  is  based  on 
solving  the  wave  equation  for  the  amplitude  of  the  laser  light  using  the  spectral  method.  This 
method  allows  us  to  model  the  interaction  of  laser  light  with  biological  cells  whose  physical 
size  can  vary  over  a  wide  range  of  spatial  length  scales,  from  sub-wavelength  dimensions  up  to 
a  few  orders  of  magnitude  larger  than  the  laser  wavelength. 


Maxwell’s  equations  in  a  nonconducting  medium  take  the  following  form  [19]  (in  Gaussian 
units): 


(1) 

(2) 

(3) 

(4) 


V  ■  e£ 
V  B 
„  -  1  dB 

c  dt 

_  juedE 
V  x5---r- 
c  dt 


=  0, 
=  0, 

=  0, 
=  0, 


where  E  and  B  are  the  electric  and  magnetic  fields,  respectively;  c  is  the  speed  of  light  in 
vacuum;  and  /i  and  £  are  the  magnetic  permeability  and  the  dielectric  constant  (or,  electric 
permittivity)  of  the  medium,  respectively.  By  taking  the  curl  of  Eq.  (3)  and  substituting  for 
V  x  B  from  Eq.  (4),  we  deduce  the  following  wave  equation  for  the  electric  field;  namely, 

(5)  ^|f-V;£  +  V(V.£)=0, 

where  we  have  used  the  vector  identity 

Vx(Vxi)  =  V(V -A)  -  V2A, 

and  the  fact  that  the  operators  Vx|  and  |Vx  are  equivalent  (i.e.,  V  and  ^  commute)  as 
space  and  time  coordinates  are  independent.  Here,  V2  =  92/9x2  +d2/dy2  +  d~/dz2  denotes  the 
Laplacian  in  a  three-dimensional  Cartesian  geometry. 


For  an  electromagnetic  wave  with  high  frequency  ff>o  say,  the  electric  field  vector  takes  the 
form 

(6)  E  =  £(x,y,z,t)exp(-/tOoO  +  c.c., 

where  c.c.  denotes  the  complex  conjugate  of  the  first  term  on  the  right-hand  side  of  this 
equation.  Here,  a>o  can  be  interpreted  as  the  carrier  frequency  of  the  electromagnetic  wave.  It 
is  assumed  that  the  amplitude  E(x,y,z,t)  varies  on  a  characteristic  time  scale  that  is  much 
larger  than  the  light  period  2jt/coo.  In  view  of  this  assumption,  enveloping  Eq.  (5)  with  respect 
to  time  (viz.,  inserting  Eq.  (6)  into  Eq.  (5)  and  ignoring  the  small-order  term  ^  ), 

reduces  the  order  of  the  time  derivative  from  two  to  one.  In  consequence,  one  obtains 


(7) 


2i^-E/u~  +  V2£  +  ^eME  -  V(V-E)  =  0. 
c2  dt  cl 
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As  a  further  simplification,  we  develop  a  wave  equation  model  for  a  one-component  field  (i.e., 
a  scalar  wave  equation  model).  For  the  one-component  field  case,  we  consider  the  scattering  of 
laser  light  by  a  biological  cell  (tissue)  in  two  and  three  spatial  dimensions.  A  2-D  scalar  wave 
equation  can  be  obtained  from  Eq.  (7)  by  assuming  that  the  biological  cell  is  inhomogeneous 
(in  electric  permittivity)  in  the  x-y  plane  only  and  that  the  electromagnetic  wave  is  ^-polarized 
[viz.,  E  —  E(x,y,t)ez  (where  ez  is  a  unit  vector  normal  to  the  x-y  plane)  so  that  the  electric  field 
polarization  vector  is  perpendicular  to  the  plane  of  inhomogeneity  of  the  biological  cell].  For 
this  case,  the  last  term  on  the  left-hand  side  of  Eq.  (7)  is  zero,  and  one  obtains  a  scalar  wave 
equation  for  the  electric  field  amplitude  ofthe  form£'(x,y,t): 

(8)  +  V2E +  ^jEjuE  =  0, 

at  c z 

where  V \  =  d2/dx2  +  d2/dy2  is  the  Laplacian  in  two  spatial  dimensions. 

The  scalar  field  approximation  can  be  also  considered  in  three  spatial  dimensions.  The  form  of 
the  wave  equation  for  the  3-D  scalar  field  approximation  is  similar  to  that  of  Eq.  (8),  except 
that  the  Laplacian  V2  has  three  spatial  components.  For  the  case  of  three  spatial  dimensions, 
the  last  term  on  the  left  hand  side  of  Eq.  (7)  may  not  be  negligible.  The  validity  of  such  an 
approximation  will  be  tested  by  comparing  simulation  results  with  exact  analytical  solutions 
(Mie  theory)  for  scattering  from  a  homogeneous  spherical  dielectric. 

For  a  biological  cell,  we  assume  fi  =  1  (non-magnetic),  and  that  the  dielectric  constant  can  be 
represented  as  a  sum  of  two  terms: 


(9)  £  =  e0  +  Ae(x,>’,z). 

Here  £o  is  the  dielectric  constant  of  the  background  fluid  (fluid  surrounding  the  biological  cell) 
that  is  considered  homogeneous,  and  (Eo  +  A e(x,y,z))  is  the  dielectric  constant  of  a  cell  that 
consists  of  cytoplasm,  nucleus,  and  other  internal  structures  and  organelles.  We  note  that  in 
Eq.  (9)  that  £  is  real  since  only  light  scattering  (and,  not  absorption)  is  considered. 
Furthermore,  measurements  of  the  optical  properties  of  Erwinia  herbicola  bacteria  [20]  have 
shown  that  the  imaginary  part  of  the  complex  index  of  refraction  N  =  \/e  is  approximately 
four  orders  of  magnitude  smaller  than  the  real  part,  implying  that  the  extinction  (or, 
absorption)  of  light  in  the  biological  cell  is  negligible  relative  to  the  light  scattering.  Hence,  an 
arbitrary  cell  model  can  be  constructed  by  simply  assigning  real  dielectric  constant  (or,  real 
electric  permittivity)  values  to  each  cell  component.  For  the  cells  considered  in  this  report,  the 
perturbations  Ae  are  assumed  to  be  small  (viz.,  |  Ae|  -C  Eq). 


It  is  readily  seen  that  the  first  term  on  the  left-hand  side  of  Eq.  (8)  accounts  for  time  variations 
of  the  electric  field  amplitude.  In  the  case  of  interest  to  us,  the  biological  cell  properties  are 
considered  to  be  stationary.  Non-stationarity  can  only  play  a  role  during  the  propagation  of  the 
electromagnetic  wave  front  through  the  biological  cell.  Hence,  £  can  be  replaced  by  £o  in  the 
time  derivative  term  in  Eq.  (8),  allowing  us  to  rewrite  this  equation  as  (for  2-D  and  3-D  cell 
geometries) 


(10) 


coq  dE 


2/^Eq^  +  V2E  +  ^e  0E  =  —  ^ f(t)AeE . 
cL  dt  cL 


m 
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The  time  function  /(/)  is  introduced  on  the  right-hand  side  of  Eq.  (10)  to  model  an  adiabatic 
“turn  on”  of  the  perturbations  with  the  property  that  f(t)  =  0  at  t  =  0  and  /(/)  — ►  1  as  t  -> 

Fig.  1  shows  the  time  function  f(t)  that  was  used  for  the  simulations  in  this  report.  This  time 
function  f(t)  allows  one  to  define  an  easy  problem  for  which  we  know  the  solution  (e.g.,  when 
f{t)  =  0  at  t  —  0  corresponding  to  the  case  of  an  electromagnetic  wave  propagating  through  a 
homogeneous  medium  with  no  scatterer),  and  a  “path”  between  this  easy  problem  and  the  hard 
problem  that  we  actually  wish  to  solve  (viz.,  electromagnetic  scattering  from  the  biological  cell 
which  is  fully  “turned  on”  in  the  simulation  when  f(t)  — ►  1  for  large  t).  Hence,  the  solution  to 
the  easy  problem  is  gradually  transformed  to  the  solution  of  the  hard  problem  by  tracing  this 
path  (which  is  defined  by  f(t)  increasing  monotonically  from  0  at  t  =  0  to  1  at  large  t). 

Though  the  unit  of  t  is  arbitrary,  for  the  simulations  conducted  in  this  report  we  will  normalize 
the  time  steps  so  that  t  is  measured  in  units  of  picoseconds  (see  details  in  the  next  section). 

The  left-hand  side  of  Eq.  (10)  accounts  for  the  field  propagation  in  a  homogeneous  medium, 
and  the  interaction  with  inhomogeneous  perturbations  (biological  cell)  is  described  by  the 
right-hand  side.  There  are  no  restrictions  in  our  model  on  the  scale  length  of  spatial  variations 
in  Ae.  However,  inhomogeneities  in  the  scattering  medium  on  a  scale  length  much  smaller 
than  the  laser  wavelength  cannot  be  properly  resolved  in  the  scattered  light  intensity  pattern. 
Therefore,  we  assume  that  the  characteristic  length  scale  of  the  spatial  inhomogeneity  for  the 
perturbation  Ae  is  larger  than  the  laser  wavelength  Xm  =  c/ (ro0  v^o)  *n  *be  background  (or, 
surrounding)  medium.  Inhomogeneities  on  a  scale  length  smaller  than  the  wavelength  X,„  (for 
example,  the  interface  between  the  cell  and  the  surrounding  fluid,  or  cell  wall)  are  treated  as 
sharp  boundaries  in  our  method.  This  approach  is  consistent  with  the  fact  that  light  scattering 
off  an  interface  much  narrower  than  the  wavelength  is  similar  to  scattering  off  a  sharp 
boundary. 

Application  of  Spectral  Method 

The  wave  equation  [Eq.  (10)]  is  solved  in  a  three-dimensional  domain  defined  as  follows: 

0  <  .v  <  Lx,  —Ly/2  <y<  Ly/2  and  -4/2  <  z  <  L-J2  (see  Fig.  2).  The  incident 
electromagnetic  wave  is  assumed  to  be  propagating  in  the  .x-direction.  Eq.  (10)  will  be  solved 
using  the  spectral  method  [18].  In  this  method,  the  electric  field  amplitude  E  is  expanded  as  a 
Fourier  series  in  the  transverse  (y  and  z)  directions  as  follows: 

(11)  e  =  Emm'  (*> exp  (imkyy + im'k*z)  > 

m  m' 

where  ky  =  2n/Ly  and  kz  =  2n/Lz  are  the  wavenumbers  in  the  y-  and  z-directions,  respectively. 
We  now  introduce  the  following  dimensionless  variables: 

T  =  Proof /eo, 

X  =  kx, 

Y  =  ky, 

(12)  Z  =  kz, 

where  k  =  2n N/Xo,  with  h  =  c/roo  being  the  laser  wavelength  in  vacuum  and  N  =  y'eo  being 
the  background  medium  refractive  index  in  the  simulation  region;  and,  p  is  a  dimensionless 
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parameter  which  is  used  to  renormalize  the  time  T  to  some  convenient  unit  (picoseconds,  for 
example).  Using  these  non-dimensional  spatial  and  temporal  variables  and  the  spectral  form 
for  the  electric  field  amplitude  E  exhibited  in  Eq.  (1 1),  Eq.  (10)  can  be  re-cast  in  the  following 
form: 

(13)  (2i^~  +  ~  +4w)  =  NL(Emm',T), 

where 

NL(Emm>,  T)  =  —Aef(t)Emm’, 

(14)  K2Xmm,  =  eo  -(mky/k)2-(m'kz/k)2. 


The  variation  of  the  electric  field  amplitude  E  in  the  longitudinal  (x)  direction  is  approximated 
using  Chebyshev  polynomials  Tn  (x)  of  the  first  kind,  which  constitute  an  orthogonal  family  on 
the  interval  - 1  <  x  <  1  with  a  weight  function  1  / \/l  —  x2.  To  facilitate  this  approximation,  the 
extent  of  the  simulation  region  in  the  direction  of  laser  light  propagation  [0,Lr]  is  mapped  into 
the  interval  [—1,1].  Now,  Emmi(X,  T)  in  Eq.  (13)  is  approximated  as 

Nx+\ 

(15)  Emm1  (x,r)  =  ^fli(7’)r,_1(z). 

k=  1 

The  use  of  Chebyshev  polynomials  evaluated  at  collocation  (grid)  points  Xj  given  by 
Xj  =  -  cos  ( (j  -  l)n/Nx),  where  j  =  1 , 2, . . .  ,NX  +  1  and  Nx  is  the  total  number  of  points  in  the 
x-direction  allows  a  fast  Fourier  transform  (FFT)  to  be  used  to  evaluate  a k(T)  in 
0((NX  +  1 )  log (Nx  +  1))  operations.  Finally,  Chebyshev  polynomials  permit  accurate 
determination  of  the  spatial  derivative  d2Emm'  /dX2  in  Eq.  (13). 

The  spectral  method  used  here  is  more  accurate  per  nodal  unknown  than  any  local  method 
(e.g.,  finite  difference,  finite  volume,  or  finite  element  method).  For  problems  with  sufficiently 
smooth  solutions  and  benign  (e.g.,  periodic)  boundary  conditions,  the  spectral  method  results 
in  a  considerable  computational  efficiency. 


Time-stepping  Procedure 

We  use  the  implicit  midpoint  rule  [21]  for  the  time-stepping  of  Eq.  (13).  The  implicit  midpoint 
method  solves  a  differential  equation  of  the  following  form 

it  =  p{u)u, 

where  p(u)  is  a  given  function  of  the  quantity  u.  For  a  step  in  time  u"  —>  w"+l , 

t„  — » tn+i  =tn  +  x,  where  x  is  the  time-step  size  and  un  is  the  quantity  u  at  t  =  tn,  the  implicit 

midpoint  rale  reads 

(16)  u"+i  =un  +  xp(u)u, 
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where  u  =  (un+i  +«")/ 2  is  the  midpoint  value  of  w.  A  solution  u"  of  Eq.  (16)  is  bounded  as 
n  — ►  oo  [22]. 


Now,  it  is  convenient  to  rewrite  Eq.  (13)  as 


(17) 


dEmm' 

dT 


i  (  32  2  ^(^m>,r)\ 

2p  v^2  w  y 


In  order  to  advance  Eq.  (17)  one  time-step,  the  forward  Euler  scheme  is  used  to  provide  a 
provisional  estimate  E*^^  for  E"+l  as 


(18) 


*{n+\) 


dEmn 

_  r?n  mn 


which  can  be  used  with  E"tm,  to  obtain  an  estimate  for  the  midpoint  value 

vn  \(r n  ,C’*(n+1))- 

^mm'  ~  2  [rnim1  )  ‘ 


(19) 


fn  ^  i  ^ 

mm'  ~  ‘-'mm'  '  qp 


a2£i 


ar2 


'-+K2Xmm,E’;wi,-NL(E"mm,,n 


where  Eq.  (17)  has  been  used  to  estimate  dEmm-/dT\n.  From  the  implicit  midpoint  rule 
[cf.  Eq.  (16)],  we  get 


(20) 


ix 


rn+l _ I _ 

c«»/  umm'  ~  2p 


32  2  _  NL(Emm'Jn) 

dX2  +  *W  £mm/ 


Substituting  Eq.  (19)  into  Eq.  (20),  assuming  that  NL(Emm'.  T")  —  NL(E"nm,,Tn),  and 
re-arranging  the  resulting  equation  we  get 


(21) 


rn+>  —  pn  j _ 

“mm'  «im  1  ^  o 


It 

4P 


a2 

dx 2 


+a:2 


Xmm' 


ESI  ' 


•,7-"), 


which  can  be  rewritten  as  follows 

(22)  (aA2  _  '£"™'  =  +NL{Emm',T  ), 

where  X,  =  —KXmm,  —  /4p /x.  The  length  of  the  simulation  region  in  the  direction  of  light 
propagation  is  mapped  into  interval  [—  1 , 1]  and,  therefore,  Eq.  (22)  is  a  Dirichlet  problem  for  a 
second-order  elliptic  operator  in  the  interval  [—1,1],  which  we  may  solve  efficiently  using  the 
Tau  method  [18]. 


Boundary  Conditions 

Among  the  issues  facing  the  implementation  of  the  numerical  algorithms  for  simulation  of 
electromagnetic  fields  in  complex  environments  is  the  proper  truncation  of  the  computational 
domain.  In  our  situation,  the  problem  to  be  simulated  corresponds  to  an  open-region  problem, 
and  the  numerical  implementation  requires  a  proper  treatment  of  the  grid  boundaries.  Due  to 
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the  properties  of  the  Fourier  transform,  the  imposed  boundary  conditions  are  implicitly 
periodic  in  the  transverse  (Y  and  Z)  directions.  The  boundary  conditions  in  Y-direction 
correspond  to  transparent  boundaries  at  X  =  X/  —  0  (left  side  domain  boundary)  and 
X  =  Xr  =  k§Lx  (right  side  domain  boundary)  for  both  the  incident  and  scattered  waves  [23]: 


BEm 


dX 


x, 


(23) 


dEn 


dX 


xr 


iKxmm'  (2 E^,(T)  —  Emm'(Xl,T)+ 

9  [E2'(E)-Emm'(XhT)]\, 


&xmm-  W 

=  iKxmm'  ( 1  +  Xt  )  (Xn  T), 

\  KXmm'  61 


where  E^m,(T)  is  the  Fourier  component  of  the  incident  electric  field  amplitude.  The 
transparent  boundary  conditions  imposed  here  in  the  longitudinal  direction  ensure  that 
spurious  reflections  from  the  grid  boundaries  at  X/  and  Xr  are  small  enough  so  that  the  solution 
is  not  contaminated. 


The  Shape  of  the  Incident  Beam 


Two  kinds  of  incident  beam  shapes  are  used  in  our  simulations.  Firstly,  the  simplest  case  is 
that  of  an  incident  plane  wave.  Secondly,  we  consider  the  Gaussian  beam  [24] 


E0{x  =  0,y,z,t)  =  E0(t)  7^7-Texp 


y2+z21 


(24) 


x  exp 


W{x ) 

■jkx-jk' 


2W2{x)\ 


2  R2{x) 


+m 


where  W (x)  and  R(x)  are  the  beam  radius  and  the  wavefront  radius,  respectively,  which  take 
the  following  forms: 


(25) 

W{x) 

=  CT 

(26) 

R(x) 

=  X 

\  x  J 

(27) 

«*) 

X 

—  arctan  — , 

Er 

where  o  is  the  minimum  radius  of  the  beam  at  the  best  focus  position  and  Lr  is  the  Rayleigh 
length,  which  is  the  distance  from  the  best  focus  position  to  where  the  beam  radius  increases 
by  a  factor  of  2.  The  minimum  radius  of  the  beam  is  related  to  the  Rayleigh  length  by 


(28) 


a  = 
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where  X  is  the  wavelength.  The  number  /  is  defined  to  characterize  the  minimum  radius  of  the 
beam  0  =  fk  at  the  best  focus  position  and  the  /-number  is  also  related  to  the  Rayleigh  length 
as  Lr  =  2n f2X.  The  Gaussian  profile  given  in  Eq.  (24)  is  a  reasonable  realization  of  an  actual 
incident  laser  beam  as  it  might  be  found  in  an  experimental  setup,  while  the  plane  wave  is  a 
idealized  profile  for  verifying  the  correctness  of  the  simulation  code,  as  it  allows  for  a  direct 
comparison  with  analytical  solutions  (e.g.,  Mie  theory). 

Far-field  Calculation 

An  appropriate  presentation  of  numerical  data  is  very  important  for  the  comparison  with 
experimental  results  and  theoretical  models.  Typical  experimental  observations  usually  only 
measure  the  characteristics  of  the  laser  light  well  outside  the  interaction  region  (the  region 
where  the  laser  light  interacts  with  the  scatterer  such  as  a  biological  cell).  Therefore,  the  laser 
field  profile  far  from  the  scattering  object  (so-called  far-field  region)  should  be  calculated  from 
the  field  distribution  in  the  interaction  region  (so-called  near-field  region)  which  can  be 
obtained  through  our  numerical  simulations.  Using  the  near-field  data  obtained  from  our 
model  simulations,  this  near-to-far-field  transformation  efficiently  and  accurately  calculates 
the  far-field  scattering  response  of  an  illuminated  bio-target  (biological  cell)  without  the  need 
to  extend  the  simulation  domain  to  the  far  field  to  obtain  the  far-field  data. 

An  asymptotic  expression  for  the  electric  field  of  the  scattered  light  at  a  far-field  point  P  (see 
Fig.  2)  is  given  by  [25]: 


m»  = 


Es(x,y,z,t)Vn 


exp  (ikp(x,y,z)  —  in/4) 


(29) 


2nkp{x,y,z) 


2nkp(x,y,z) 
exp(ikp(x.\\z)  -  in/4)'V„Es(x.y.z,t ) 


where  Es  is  the  electric  field  at  the  emitting  surface  S,  E(P,t)  is  the  electric  field  at  a  far-field 
point  P,  V„  denotes  the  derivative  along  the  outwardly-directed  normal  vector  to  the  emitting 
surface,  and  p  is  the  vector  from  the  point  at  the  emitting  surface  with  coordinates  (x,y,z)  to 
the  far-field  point  with  coordinates  (x0,yoFo)-  Furthermore,  let  p0  be  the  postion  vector  from 
the  origin  of  the  coordinate  system  to  the  far  field  point  [viz.,  po  =  (xo,yo,zo)]-  In  order  to 
evaluate  the  surface  integral  of  Eq.  (29),  we  use  a  Cartesian  geometry  and  let 


u  =  (l/2nkp(x,y,z))exp(ikp(x,y,z)  -  in/ 4). 


Then  the  surface  integral  of  Eq.  (29)  takes  the  form 

(30)  E(P,t)  =  ^dS-{EV„u-uV„E). 

From  Fig.  2,  one  can  see  that  the  planes  ABCD  and  EFGH  coincide  with  the  inlet  (x  =  0)  and 
outlet  (x  =  Lx)  boundaries.  For  forward  scattering  calculations,  we  consider  only  the 
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contributions  due  to  plane  EFGH,  while  for  backward  scattering  calculations,  we  consider 
only  the  contributions  due  to  plane  ABCD.  Hence,  the  integral  of  Eq.  (30)  reduces  to 


(31) 


Em=hjldyd2(ETx-ui*) 


The  above  result  can  be  expressed  in  a  somewhat  more  informative  form  if  we  rewrite  p  in 
terms  of  spherical  coordinates  (p,0,<j)),  where  p  =  |p|.  Also,  we  let  po  =  |po|  denote  the 
distance  from  the  origin  of  the  coordinate  system  to  the  far  field  point.  Using  the  following 
mathematical  relationships 


(32) 
we  have 


(33) 


3p/dx 

du/dx 

P 

lim  u 

p — too 


du 
lim  — 

p-*~  dp 


=  -x/p  —  —  cos  9, 

=  (3zz/3p)(3p/3x)  —  —  cos9(3w/3p), 

~  Po -xcos0  —  ysin0cos<j>-zsin0sin<|>, 

=  - — - — exp(z'£po  —  z'kxcos0  —  z'Aysin0cos<{> 
Inkpo 

—ikz  sin  0  sin  (j)  —  z  Jt/4) , 
ik 

=  — -—exp  (ikpo  —  ikx  cos  9  —  iky  sin  0  cost)) 
2jtA:po 

— /fe  sin  0  sin  0  —  m/4) , 


E(P,t)  = 


1  1 


;// 


2 i 2nkpo 
x  exp  (—z'fct  cos  0 


dydz 


•zzr  n  dE 
—ikE  cos  9  —  — 

OX 


■  iky  sin  0  cos  <|)  —  ikz  sin  0  sin  <|>) . 


By  introducing  kx  =  £cos0,  ky  =  &sin0cos<|),  kz  =  £sin0sin<|>,  and  taking  into  account  that,  in 
the  case  of  forward  scattering,  E  ~  exp  ( ikxx ),  one  obtains  from  Eq.  (33)  that 


(34) 


E{P,t)  =  CXP ^fP,+  [~2ikxexp(-ikxx)]E(Lx,kv,k:). 

Qiikpo 


Similarly,  in  the  case  of  backward  scattering,  we  have 
(35) 


E{P,t)  =  CXP in^4‘>  [2ikxexp(-ikxx)}E(0,ky,kz). 


4nkpo 


Eq.  (34)  for  E(P,t)  has  been  used  in  our  numerical  simulations  to  calculate  the  forward 
scattered  light  intensity  in  the  far-field  region.  The  forward  emitted  power  at  a  far-field  point  is 
then  given  by 


(36) 


Q=^Po\E(P,t)\2  =  \ikxE(Lx, ky, kz) |2 . 
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Code  Implementation  and  Verification  _ 

In  the  previous  section,  we  presented  the  theoretical  (albeit  approximate)  model  equations 
upon  which  our  numerical  simulations  are  based.  More  specifically,  we  derived  an 
approximate  scalar  3-D  wave  equation  from  Maxwell’s  equations,  described  the  numerical 
algorithm  which  we  will  use  to  solve  this  equation,  and  showed  how  to  compute  the  scattering 
pattern  in  the  far-field  region  using  a  near-field  to  far-field  integral  transformation. 

For  the  model  equations  described  in  the  previous  section  for  both  2-D  and  3-D  cell 
geometries,  a  computer  code  has  been  implemented  to  solve  these  model  equations  using  the 
spectral  method.  Since  the  code  is  new,  it  needs  to  be  verified  before  being  used  for  any 
applications.  In  this  section,  results  for  both  2-D  and  3-D  cell  geometries  will  be  presented  and 
verified  against  predictions  obtained  from  linear  perturbation  theory  and  Mie  theory, 
respectively. 


Linear  Perturbation  Theory  Applied  to  2-D  Scattering 


In  this  subsection,  an  analytical  stationary  solution  to  Eq.  (10)  is  obtained  using  the  linear 
perturbation  method.  In  the  stationary  regime,  the  time  derivative  term  in  Eq.  (10)  is 
identically  zero.  Now,  consider  the  scalar  wave  equation  [Eq.  (10)]  in  two  spatial  dimensions. 
By  assuming  small-angle  scattering,  the  field  amplitude  E  can  be  written  as 
E  =  Ea(x,y)exp(ikmx),  where  k„,  =  2n/Xm  =  (coo/c)\/^o  is  the  laser  wavenumber  in  the 
background  medium.  The  characteristic  inhomogeneity  length  scale  of  the  field  amplitude 
Ea(x,y)  in  the  direction  of  light  propagation  (x)  is  assumed  to  be  much  larger  than  the  laser 
wavelength  in  the  background  medium,  so  X.OT  |91n(E-0)/3x|“1 .  This  is  precisely  the 
condition  that  is  required  for  validity  of  the  paraxial  approximation.  Expanding  Ea  =  Eq  +  8E  , 
where  Eq  is  the  incident  light  amplitude  and  5 E  is  the  scattered  light  amplitude  (]8£|  <§C  |£o|)> 
Eq.  (10)  to  first  order  in  5 E  and  Ae  assumes  the  form 


(37) 


„  ,  95 E  92S E 

2lkm  dx  +  dv2 


cd; 


—j-AeEq  =  0. 


After  performing  a  Fourier  transformation  in  the  transverse  (y)  direction,  we  obtain  the 
following  differential  equation  for  the  Fourier  component  of  the  electric  field  8 E(x.ky)  : 


(38)  2  ikm  d?>E^  —  -  E  {x,  ky)  +  ^  Ae(x,  ky)E0  =  0, 

where  A e(x,kv)  is  the  Fourier  transform  of  the  permittivity  perturbation  Ae(x,y)  in  the 
y-direction.  The  solution  to  Eq.  (38)  takes  the  form  of  the  following  integral  transform 

(39)  5£M,)  =  Ae(,',*,)exp  ■ 

The  constant  of  integration  in  Eq.  (39)  is  chosen  so  that  the  field  perturbation  8£  is  zero  at  the 
boundary  x  =  0. 
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For  simplicity,  consider  a  square  perturbation  of  the  dielectric  constant  A e(x,y)  equal  to  Aeo  in 
the  region  -Ro  <  x -xo  <  Ro  and  -Ro  <y<  Ro,  and  zero  elsewhere.  By  using  the  assumption 
that  the  incident  light  has  a  very  small  divergence  angle  ( Eo(Lx,ky )  ~  0  for  ky  ^  0),  we  may 
combine  Eqs.  (39)  and  (36)  to  obtain  the  transmitted  power  per  unit  angle  as 


(40) 


Q  = 


8 El 
n3Lj 


1-^ 


sm 


sin 


where  k *  =  ky/ko. 

Similarly,  placing  a  square  nucleus  measuring  2Rn  by  2R„  and  dielectric  constant  perturbation 
Ae„  in  the  center  of  a  cell  with  dimensions  2Rq  by  2Ro  and  dielectric  constant  Aeo,  the 
transmitted  power  per  unit  angle  is 


As  L  — >  0,  Eq.  (41)  takes  the  form 


»  a?  (?)'©' 

The  scattering  from  these  square  perturbations  can  be  used  to  verify  the  numerical  results 
obtained  from  our  computer  code  implementing  the  spectral  method.  A  comparison  between 
our  numerical  results  and  the  theoretical  predictions  of  Eqs.  (40)  and  (41)  is  shown  in  Fig.  3. 
In  these  simulations,  the  dielectric  constant  of  the  medium  is  £o  =  1 .8225,  while  the  dielectric 
constant  perturbations  for  the  cytoplasm  and  the  nucleus  are  Aeo  =  0.0544  and  Ae„  =  0.1096, 
respectively.  The  accuracy  of  our  numerical  method  in  this  case  is  clearly  demonstrated  in 
Fig.  3,  where  the  numerical  simulation  results  show  good  agreement  with  the  theoretical 
calculations. 


Scattering  from  a  3-D  Homogeneous  Sphere 

The  diffraction  of  a  plane  monochromatic  wave  by  a  homogeneous  sphere  with  a  zero 
thickness  interface  (no  coating  on  sphere)  can  be  calculated  within  the  framework  of  Mie 
theory  [25,  26],  Mie  theory  exactly  solves  Maxwell’s  equations  in  three  spatial  dimensions  for 
all  components  of  the  electromagnetic  field  both  inside  and  outside  the  sphere.  In  Mie  theory, 
the  electromagnetic  fields  are  expressed  in  terms  of  infinite  series  expansions.  These 
expansions  can  be  significantly  simplified  when  the  distance  r  from  the  scattering  sphere  is 
much  larger  than  the  laser  wavelength  Xm  in  the  medium  surrounding  the  sphere.  When  the 
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incident  light  is  linearly  polarized,  the  scattered  light  intensity  |S|2  can  be  determined  from 

[26]  o  0.1 

E2  Yl  *T"  1  r  ,  1 

/I  v  ' 

where  7t„  =  -(l/sin0yP„/</0,  x„  =  - d2Pn/dQ 2,  and  /»„  are  the  Legendre  polynomials  of 
order  n.  The  coefficients  a„  and  b„  are  given  by 

_  wvi/^iiocXCy)  -\|/n(.t)VH(wtr) 

a"~  -c„(xX(mx)  ’ 

\|/„(mx)\j/;,(x)  -  my„(x)V„H 

(44)  \l/„(/nx)^(x)-/nc„(x)<(/wx)  ’ 

where  \|/„  and  are  the  Ricatti-Bessel  functions,  x  =  kma,  m  =  ks/km,  a  is  the  radius  of  the 
scattering  sphere,  and  km  and  ks  are  the  laser  wavenumbers  outside  and  inside  the  scattering 
sphere,  respectively. 

Figs.  4  and  5  compare  3-D  simulation  results  and  Mie  theory  predictions  for  the  scattering  of  a 
plane  wave  by  a  homogeneous  sphere  in  two  different  cases.  The  agreement  is  very  good  for 
scattering  angles  up  to  about  30  degrees.  This  good  agreement  between  our  numerical 
simulations  and  Mie  theory  is  achieved  in  spite  of  the  fact  that  the  numerical  model  solves  the 
wave  equation  in  the  scalar  approximation  only,  whereas  Mie  theory  solves  the  full  Maxwell  s 
equations  (for  all  six  components  of  the  electromagnetic  field).  We  may  conclude  that  the 
scalar  wave  equation  is  an  effective  approximation  to  the  full  Maxwell  equations  for  moderate 
scattering  angles  that  are  of  practical  importance  for  applications  such  as  flow  cytometry.  The 
scalar  wave  approximation  is  valid  because  in  the  case  of  small-angle  scattering,  the 
longitudinal  component  of  the  electric  field  E  and  V  •  E  are  both  small,  so  the  last  term  on  the 
left-hand  side  of  the  full  wave  equation  [Eq.  (7)]  can  be  neglected.  For  larger  scattering  angles, 
the  scalar  wave  approximation  breaks  down,  since  the  term  with  V  •  E  becomes  comparable  in 
magnitude  to  the  second  and  third  terms  on  the  left-hand  side  of  Eq.  (7). 


Applications  and  Discussion _ _ _ _ 

The  previous  two  sections  dealt  with  the  development  of  an  algorithm  upon  which  our 
simulations  are  based  and  the  implementation  and  verification  of  the  computer  code.  In  this 
section,  we  present  some  preliminary  two-dimensional  and  three-dimensional  simulation 
results  of  laser  light  scattering  from  biological  cells  to  demonstrate  some  of  the  potential 
applications. 

The  biological  cell  with  typical  size  in  the  range  of  10-30  /xm  consists  of  various  organelles 
immersed  in  a  cytoplasm,  all  of  which  are  enclosed  within  a  cell  wall  [27].  The  largest 
organelle  in  the  cell  is  usually  the  nucleus,  whose  size  ranges  typically  between  3  and  10/an 
[28].  Since  we  intend  to  model  laser  light  scattering  in  biological  cells,  we  require  values  for 
the  index  of  refraction  of  each  cell  component.  The  cell  is  defined  by  assigning  electric 
permittivity  values  (or,  equivalently,  refractive  index  values)  to  each  cell  component.  The 
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values  listed  in  Table  1  were  taken  from  the  literature,  and  are  used  in  all  simulations  in  this 
report. 


Cell  component 

Refractive  index 

Reference 

Surrounding  medium 

1.35 

[29] 

Cytoplasm 

1.37 

[30] 

Nucleus 

1.39 

[30] 

Mitochondria 

1.42 

[31] 

Table  1:  The  values  for  index  of  refraction  of  various  cell  components  which  are  used  in  the  simulations.  The 
medium  surrounding  the  cell  is  assumed  to  be  water. 


Two-Dimensional  Simulation  Results 

Although  they  cannot  describe  the  “real  world”,  results  from  simulations  in  two  dimensions 
are  nevertheless  useful  in  obtaining  qualitative  predictions  for  the  behaviour  of  certain 
systems.  The  2-D  code  is  computationally  more  efficient  than  the  3-D  code,  and  is 
advantageous  in  that  memory  requirements  are  modest.  Nevertheless,  it  is  expected  that  the 
2-D  scattering  solutions  from  heterogeneous  biological  cells  will  exhibit  qualitatively  similar 
trends  to  the  more  complex  3-D  solutions.  In  the  2-D  case,  we  will  focus  on  plane-wave 
scattering  from  both  single  and  multiple  biological  cells  in  the  sampling  volume. 

In  light  scattering  from  cells,  it  is  important  to  understand  the  role  of  a  cell’s  nucleus  on  the 
scattering  intensity  pattern.  To  this  end,  the  first  set  of  simulations  was  designed  to 
demonstrate  how  the  size  of  the  nucleus  in  a  cell  contributes  to  scattering.  Three  types  of 
biological  cell  morphologies  were  chosen  for  our  2-D  simulations:  namely,  a  cell  without  a 
nucleus,  a  cell  with  a  normal-size  nucleus  of  diameter  Dn  =  3  /mn,  and  a  cell  with  a  large 
nucleus  of  diameter  D„  —  6  ym  (for  example,  it  is  known  that  cancerous  cells  are  characterized 
by  a  large  nuclear-to-cytoplasmic  volume  ratio  [32]).  For  our  simulations,  the  biological  cell 
diameter  in  all  three  cases  is  taken  to  be  Do  —  1 1  /mi,  which  is  within  the  range  of  typical 
biological  cell  sizes  cited  above  (more,  specifically,  this  size  is  typical  for  an  animal  cell,  but  is 
generally  larger  than  a  bacterial  cell).  The  dielectric  constant  (or,  electric  permittivity)  £  is 
related  to  the  index  of  refraction  N  by  e  =  N2.  Therefore,  using  the  values  from  Table  1,  we 
know  that  the  dielectric  constant  of  the  surrounding  medium  (water)  is  £o  =  1 .8225.  Similarly, 
we  find  that  for  the  cytoplasm,  £  =  1 .8769,  while  for  the  nucleus  £  =  1 .9321 .  The  incident 
laser  light  is  assumed  to  be  a  monochromatic  plane  wave  with  wavelength  Xg  =  1  jum  (near 
infrared  region).  The  dependence  of  the  scattered  light  intensity  on  the  scattering  angle  is 
shown  on  a  linear  scale  in  Fig.  6(a)  and  on  a  logarithmic  scale  in  Fig.  6(b).  The  transmitted 
component,  which  is  the  large  peak  at  0  degree  in  Fig.  6(b),  is  removed  by  setting  the  intensity 
at  0  degree  equal  to  its  neighbouring  point  intensity  in  Fig.  6(a). 

From  Fig.  6,  it  is  clear  that  the  dominant  feature  of  the  scattering  pattern  is  contributed  by  the 
forward  (small-angle)  scattering.  Three  different  regions  can  be  identified  in  the  scattering 
pattern  in  Fig.  6.  First,  the  angular  interval  from  0  to  0.25  degrees  corresponds  to  the  angular 
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width  of  the  incident  light.  The  large  scattering  intensity  over  this  range  simply  represents  the 
signal  from  the  incident  light  that  is  not  scattered.  The  second  region,  which  extends  from  0.25 
to  0o  ~  4  degrees,  contains  the  main  part  of  the  scattered  light  intensity  from  the  biological 
cell.  The  width  of  this  interval,  characterized  by  the  value  of  0o,  is  determined  by  the  ratio  of 
the  typical  cell  size  to  the  laser  wavelength  in  the  surrounding  medium.  A  cell  with  a  large 
nucleus  scatters  more  light  than  a  cell  with  a  normal  nucleus,  while  a  cell  with  a  normal 
nucleus  scatters  more  light  than  a  cell  without  a  nucleus.  This  is  consistent  with  the  predictions 
of  linear  perturbation  theory  [Eqs.  (41)  and  (42)]  and  with  the  FDTD  numerical  results  of  [12]. 
Finally,  for  angles  larger  than  90  se  4  degrees,  the  scattering  amplitude  decreases  rapidly  with 
increasing  scattering  angle,  so  that  the  large-angle  scattering  intensity  is  small  in  magnitude. 

In  many  realistic  situations,  the  laser  light  propagates  through  a  medium  (sampling  volume) 
with  many  scattering  objects  (scatterers).  For  example,  a  sampling  volume  probed  by  the  laser 
light  may  contain  many  biological  cells.  Therefore,  the  problem  of  light  scattering  by  an 
ensemble  of  identical  scatterers  is  of  fundamental  importance.  Next,  the  scattering  pattern  of  a 
scattering  volume  containing  many  scattering  objects  will  be  described  and  compared  to  the 
scattering  pattern  of  a  single  scatterer.  In  our  2-D  simulations,  the  number  of  randomly 
distributed  scattering  objects  in  the  simulation  region  was  varied  from  1  to  90.  The  scattering 
objects  are  assumed  to  be  homogeneous  spheres  each  with  a  diameter  of  A)  =  7  /mi  and 
permittivity  of  e  =  1.8769.  The  permittivity  of  the  surrounding  medium  (water)  is  assumed  to 
be  e0  =  1 .8225.  The  incident  light  is  taken  to  be  a  monochromatic  plane  wave  with  a 
wavelength  Xq  =  1  fan. 

Fig.  7  illustrates  the  typical  form  of  the  relationship  between  the  scattering  intensity  and  the 
scattering  angle  for  the  case  of  10  random  scatterers  and  of  a  single  scatterer  in  the  simulation 
domain.  It  is  readily  seen  that  the  scattering  intensity  from  10  scatterers  does  not  have  the 
regular  pattern  of  minima  and  maxima  that  is  observed  for  laser  light  scattering  from  only  a 
single  scattering  center.  Not  surprisingly,  the  scattering  intensity,  averaged  over  an  angular 
interval  of  a  few  degrees,  is  several  times  larger  for  the  10  scatterer  case  than  for  the  single 
scatterer  case.  The  increase  in  the  scattering  intensity  with  an  increasing  number  of  scatterers 
in  the  sampling  volume  is  observed  at  both  small  and  large  scattering  angles.  One  can 
conclude  from  Fig.  7  that  the  increase  in  the  number  of  scatterers  changes  the  magnitude  of  the 
scattering  intensity,  but  although  the  broad  characteristic  angular  features  of  the  scattering 
pattern  are  roughly  unchanged  compared  to  that  of  a  single  scatterer  it  is  seen  that  the  random 
distribution  of  the  multiple  scatterers  results  also  in  a  complex  microstructure  being 
superimposed  on  the  general  scattering  pattern. 

The  dependence  of  scattering  power  and  scattering  angle  on  the  number  of  scatterers  is 
illustrated  in  Figs.  8(top)  and  (bottom),  which  combine  results  corresponding  to  two  different 
lengths  of  the  simulation  region;  namely,  Lx  =  160  //m  and  Lx  =  320  fim.  Fig.  8(top)  shows  the 
fraction  of  the  scattered  power  (i.e.,  the  power  of  light  propagating  outside  the  incident  beam 
optics)  as  a  function  of  the  number  of  scatterers.  The  scattered  power  varies  from  about  4%  for 
a  single  scatterer  to  over  90%  for  90  scattering  centers.  From  Fig.  8(top),  it  is  evident  that  the 
scattered  power  does  not  depend  on  the  length  of  the  simulation  region  (sampling  volume),  but 
only  on  the  number  of  scatterers  contained  in  the  sampling  volume.  The  scattered  power  is 
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nearly  linearly  dependent  on  the  number  of  scatterers  in  the  sampling  volume  [cf.  Fig.  8(top)]. 


The  angular  width  of  the  scattered  light  can  be  characterized  by  the  average  scattering  angle 
(0),  defined  as 


(45) 


fe2Q(e)de\ 
fQ(8)de  ) 


1/2 


where  Q  is  the  transmitted  power  per  unit  angle.  Fig.  8(bottom)  shows  the  dependence  of  the 
average  scattering  angle  (9)  on  the  number  of  scatterers.  It  is  evident  from  this  figure  that  (0) 
changes  only  slightly  (by  no  more  than  about  25%)  when  the  number  of  scatterers  changes 
from  1  to  90.  The  small  dependence  of  the  characteristic  scattering  angle  (0)  on  number  of 
scatterers  is  consistent  with  the  angular  spectrum  of  the  scattering  intensity  that  showed  an 
increase  in  the  amount  of  scattering  without  essential  changes  in  the  angular  spectrum 
[cf.  Fig.  7]. 


It  has  been  demonstrated  that  the  spectral  method  can  be  applied  to  model  laser  light  scattering 
from  scattering  volumes  that  contain  a  single  biological  cell  or  a  large  number  of  biological 
cells.  In  the  case  of  laser  light  scattering  from  a  single  biological  cell,  the  most  significant 
features  of  the  angular  distribution  of  the  scattered  light  from  various  cells  occur  at  small 
angles.  Cells  containing  large  nuclei  have  more  scattering  power  at  small  angles  compared  to 
cells  containing  normal-sized  nuclei.  The  amplitude  of  the  scattered  light  decreases  even 
further  for  cells  without  a  nucleus.  In  the  case  of  scattering  from  multiple  cells,  the  simulation 
results  show  that  the  fraction  of  the  scattered  power  depends  on  the  number  of  scatterers  and 
not  on  the  size  of  the  simulation  region.  In  fact,  the  fraction  of  the  scattered  power  is  almost 
linearly  proportional  to  the  number  of  scatterers.  Finally,  the  scattering  pattern  from  many 
cells  still  roughly  retains  the  angular  distribution  characteristics  of  single  cell  forward 
scattering,  but  with  a  high  frequency  modulation  due  to  interference  effects. 


Three-Dimensional  Simulation  Results 

Since  the  three-dimensional  model  more  closely  describes  real-world  applications,  the  results 
from  the  3-D  simulations  are  very  useful  in  obtaining  quantitative  predictions  for  the  behaviour 
of  laser-cell  interactions.  In  consequence,  these  results  have  many  potential  applications  not 
the  least  of  which  is  to  allow  one  to  assess  (diagnose)  biological  cell  morphology  from  light 
scattering  measurements  such  as  those  obtained  in  flow  cytometry  measurements.  The  latter 
application  can  lead  potentially  to  improved  detection  performance  for  BW  agents. 

In  the  3-D  case,  the  incident  laser  light  profile  can  either  take  the  form  of  a  plane  wave  or, 
more  realistically,  that  of  a  Gaussian  beam.  In  this  subsection,  the  influence  of  cell  size  and 
nucleus  size  on  the  scattering  intensity  are  studied.  Figs.  9-12  show  the  angular  spectra  of 
light  scattered  by  a  concentric  spherical  cell  in  3-D  numerical  simulations.  The  angular  range 
shown  in  Figs.  9-12  is  up  to  45  degrees  away  from  the  direction  of  incident  light  propagation. 
Fig.  9  shows  the  angular  spectra  of  scattered  light  for  the  case  when  the  incident  laser  light  is  a 
monochromatic  plane  wave  with  wavelength  Xq  =  1  /un.  The  length  of  the  simulation  region  is 
40Xo  in  the  direction  of  light  propagation  (^-direction)  and  80^o  in  each  of  the  two  transverse 
directions  (y -  and  z-directions).  In  Fig.  9,  three  types  of  cells  are  used  in  our  3-D  simulations; 
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namely,  a  cell  without  a  nucleus,  a  cell  with  a  normal-sized  nucleus  (the  diameter  of  the 
nucleus  Dn  =  3  /an),  and  a  cell  with  a  large  nucleus  (the  diameter  of  the  nucleus  D„  =  6  ^m). 

The  cell  diameter  in  all  cases  is  D0  =  1 1  /«n.  The  refractive  index  of  the  surrounding  medium 
is  assumed  to  be  No  =  \/£o  —  1-35.  The  heterogeneous  cell  is  modelled  as  a  concentric 
spherical  geometry  (cytoplasm  plus  nucleus)  with  refractive  indices  of  .V  =  1 .37  and  1.39  for 
the  cytoplasm  and  nucleus,  respectively.  A  spherical  nucleus  is  located  at  the  center  of  the 
spherical  cell. 

Fig.  10  shows  the  angular  scattering  spectra  for  cells  without  a  nucleus,  but  with  different 
diameters  (sizes);  namely,  for  Z)o  =  11  A)  =  7  fira,  and  Do  —  5  fim.  The  scattering 
intensity  patterns  in  Figs.  9  and  10  have  sharp  maxima  at  zero  angle,  corresponding  to  the 
unscattered  (transmitted)  light.  In  our  simulations  with  a  plane  wave,  only  a  small  fraction  of 
light  is  scattered  because  the  cell  volume  occupies  only  a  small  fraction  of  the  simulation 
region.  Figs.  9  and  10  clearly  show  that  the  scattered  light  spectra  have  a  characteristic  width 
of  a  few  degrees  and  that  this  width  is  consistent  with  the  cell  size;  viz.,  the  smaller  the  cell 
size,  the  larger  the  angular  width  of  the  main  lobe  of  the  scattered  intensity  pattern.  In  fact,  the 
first  minimum  in  the  scattering  pattern  is  approximately  determined  by  the  ratio  of  the  incident 
laser  wavelength  and  the  diameter  of  the  cell.  For  example,  the  three  cases  of  cell  diameter 
shown  in  Fig.  10  have  minima  at  (1/11)  x  57.3°  «  5.2°,  (1/7)  x  57.3'  «  8.2°,  and 
(1/5)  x  57.3°  ~  1 1.4°,  respectively,  for  biological  cells  with  diameters  of  11,  7,  and  5 /on. 

This  observation  is  consistent  (approximately  or  better)  with  diffraction  theory  (e.g.,  the 
intensity  distribution  of  diffraction  from  a  spherical  particle  of  diameter  d  will  have  have  its 
first  minimum  roughly  at  angle  0  ~  "k/d,  where  X  is  the  wavelength  of  the  incident  light  beam 
which  is  consistent  with  the  results  shown  in  Fig.  10). 

Figs.  1 1  and  12  show  the  angular  spectra  of  scattered  light  in  the  case  when  the  incident  laser 
light  is  a  Gaussian  beam  with  a  wavelength  of  Xo  =  1.0  /im  and  optical  /-numbers  of//3 
(Fig.  1 1)  and  // 10  (Fig.  12).  The  results  of  Fig.  1 1  are  for  scattering  from  homogeneous  and 
nucleated  spherical  cells  each  with  a  diameter  of  1 1  n m,  whereas  those  of  Fig.  1 2  are  for 
homogeneous  spherical  cells  of  different  diameters.  By  comparing  Figs.  1 1  and  12  with  Figs.  9 
and  10,  it  can  be  seen  that  in  the  case  of  a  Gaussian  beam,  the  amount  of  scattered  light  is 
much  larger  than  in  the  case  of  a  plane  wave.  This  is  partially  due  to  a  stronger  concentration 
of  the  total  light  beam  on  a  cell.  The  amount  of  scattered  light  can  provide  information  about 
the  permittivity  of  a  cell.  This  information,  combined  with  the  information  on  the  cell  size 
from  the  width  of  the  angular  spectrum,  allows  a  characterization  of  the  basic  cellular  features. 

Figs.  1 1  and  12  also  show  that  the  transmitted  component  of  the  incident  light  dominates  the 

angular  spectra  of  the  scattered  light.  In  order  to  analyze  properly  the  scattered  light  intensity 

pattern,  the  transmitted  component  of  incident  light  must  be  removed.  If  the  incident  laser 

light  impinges  directly  on  a  biological  cell  and  a  detector  is  placed  in  the  forward  direction  to 

measure  the  forward  scattering  light  intensity  pattern  (i.e.,  the  scattering  intensity  at  small  « 

scattering  angles),  then  it  seems  that  a  straightforward  solution  to  the  problem  would  be  to 

simply  insert  an  small  opaque  disc  in  front  of  the  detector  to  block  the  transmitted  component 

of  the  incident  light.  However,  the  problem  with  this  simple  solution  is  that  diffraction  effects 

arising  from  the  disc  edge  will  distort  the  scattering  intensity  pattern  in  the  forward  direction. 
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A  much  more  elegant  solution  to  this  problem  would  be  to  use  a  metal  clad  leaky  waveguide  to 
detect  the  scattered  intensity  spectrum  as  described  recently  in  [33].  The  basic  idea  here  is  to 
measure  the  light  scattering  pattern  from  a  bio-target  resulting  not  from  the  incident  laser  light 
illumination,  but  rather  from  an  evanescent  wave  illumination  generated  as  a  guided  optical 
wave  in  a  metal  clad  waveguide  device. 

Summary  and  Conclusions _ _ 


This  report  presents  some  preliminary  results  from  numerical  and  theoretical  studies  of  the 
effects  of  permittivity  heterogeneities  on  light  scattering  properties  from  biological  cells  and, 
in  particular,  of  how  two  major  parts  of  the  cell  contribute  to  this  scattering.  The  spectral 
method  is  applied  to  solve  two-  and  three-dimensional  scalar  wave  equations  that  are  used  to 
approximate  the  full  set  of  Maxwell’s  equations.  We  have  implemented  both  2-D  and  3-D 
solvers  for  these  scalar  wave  equations  and  applied  them  to  the  prediction  of  the  scattered 
intensity  patterns  from  biological  cells  containing  heterogeneous  organelles. 

The  accuracy  of  the  2-D  scalar  solver  has  been  verified  by  comparison  with  linear  perturbation 
theory.  Furthermore,  Mie  theory  has  been  applied  to  verify  the  3-D  scalar  solver  for  a 
homogeneous  spherical  scatterer.  A  comparison  of  the  scattered  intensity  for  a  homogeneous 
spherical  scatterer  also  demonstrated  that  the  3-D  scalar  wave  equation  provides  a  good 
approximation  to  the  full  set  of  Maxwell’s  equations  for  light  scattering  at  small  to  moderate 
angles  (up  to  about  35°). 

The  simulation  code  can  be  used  to  model  light  scattering  from  both  a  single  cell  or  a  large 
number  of  cells  in  the  sampling  volume  (simulation  region).  The  spectral  method  used  to  solve 
the  scalar  wave  equations  here  allows  a  very  high  spatial  resolution,  limited  only  by  the 
dimensions  of  the  simulation  region  (box)  in  the  transverse  direction.  For  example,  if  the  laser 
wavelength  Xo  is  1.0  fim  and  Ly  =  LZ  —  80Xo,  then  the  resolution  is  «  0.53  degrees,  which  is 
comparable  to  the  resolution  attainable  in  experimental  measurements  using  a  flow  cytometer. 
The  methodology  accurately  models  cell  inhomogeneities  on  spatial  scales  larger  than  the 
laser  wavelength,  and  treats  inhomogeneities  with  a  spatial  scale  much  smaller  than  the  laser 
wavelength  as  sharp  boundaries.  The  computational  requirements  for  application  of  a  spectral 
method  (using  the  fast  Fourier  transform  algorithm)  to  the  solution  of  a  scalar  wave  equation 
are  much  lower  than  those  for  solution  of  the  full  set  of  Maxwell’s  equations  using  a  numerical 
method  such  a  finite-difference  time-domain  (FDTD)  method  which  would  involve  discretizing 
Maxwell’s  curl  equations  in  space  and  time  to  give  a  set  of  explicit  finite-difference  equations 
that  would  need  to  be  stepped  in  time  (with  the  electric  and  magnetic  fields  at  each  grid  point 
in  the  simulation  domain  being  alternately  updated).  However,  the  price  paid  for  this 
computational  efficiency  is  that  the  current  model  only  solves  a  scalar  wave  equation  so  all 
polarization  information  is  lost.  Furthermore,  the  method  cannot  resolve  explicitly  structures 
(e.g.,  cell  membrane)  on  scales  smaller  than  the  probing  laser  wavelength. 
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“Start  up"  Profile 


Figure  1:  The  profile  of  time  function  f(t)  which  is  used  to  “turn  on"  the  perturbations.  Here  the  simulation  time  is 
normalized  and  a  stationary  solution  is  generated  in  the  time  interval  from  40  to  50  ps. 
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Figure  2:  The  geometry  of  the  simulation  region. 
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Angle[degrees] 

Figure  3:  A  comparison  of  the  scattering  intensity  obtained  from  both  numerical  simulations  and  analytic  solutions 
[Eq.  (42)]  for  three  types  of  square  cells:  for  a  cell  without  a  nucleus— the  dash-dotted  line  (simulations)  and  o 
(analytical):  for  a  cell  with  a  normal  nucleus— the  solid  line  (simulations)  and  4-  (analytical):  and  for  a  cell  with  a  large 
nucleus — the  dashed  line  (simulations)  and  A  (analytical).  The  diameter  of  the  cell  isDo  =  11  pm,  the  width  of  the 
normal  nucleus  is  D„  =  3  pm,  and  the  width  of  the  large  nucleus  is  Dn  =  6  pm.  The  dielectric  constant  perturbation 
for  the  cytoplasm  is  Aeo  =  0.0544,  while  that  for  the  nucleus  is  Ae„  =  0. 1096. 
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scattered  intensity  [a.u] 


Figure  4:  A  comparison  of  the  intensity  of  scattered  light  from  a  homogeneous  sphere  in  three  spatial  dimensions 
obtained  from  simulations  (solid  line)  and  from  Mie  theory  (dotted  line).  The  diameter  of  the  sphere  is  Dq  =  11  pm. 
The  dielectric  constant  outside  the  sphere  is  eo  =  l  and  inside  the  sphere  ise=  1 .0609.  In  the  simulation  result,  the 
transmitted  component  (at  (f)  has  been  removed  by  setting  the  intensity  here  equal  to  its  neighbouring  point. 
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Figure  5:  Another  comparison  of  the  intensity  of  scattered  light  from  a  homogeneous  sphere  in  three  spatial 
dimensions  obtained  from  simulations  (solid  line)  and  from  Mie  theory  (dotted  line).  The  diameter  of  the  sphere  is 
Do  =  11  pm.  The  dielectric  constant  outside  the  sphere  is  Eo  =  1.8225  and  inside  the  sphere  is  e  =  1.8769.  In  the 
simulation  result,  the  transmitted  component  (at  CP)  has  been  removed  by  setting  the  intensity  here  equal  to  its 
neighbouring  point. 
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Figure  6:  Scattering  intensity  from  a  two-dimensional  single  cell  with  diameter  Do  =  11  pm  on  a  linear  scale  (top) 
and  a  semi-logarithmic  scale  (bottom).  The  dash-dotted  line  corresponds  to  a  cell  without  nucleus,  the  solid  line 
corresponds  to  a  cell  with  a  normal  nucleus  of  diameter  D„  =  3  pm,  and  the  dashed  line  corresponds  to  a  cell  with  a 
large  nucleus  of  diameter  Dn  =  6  pm.  The  transmitted  component,  which  is  the  large  peak  at  CP  (bottom),  is 
removed  by  setting  the  intensity  here  equal  to  a  neighbouring  point  intensity  (top). 
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as  a  function  of  scattering  angle  for  10  scatterers  (dotted  line)  and  for  a  single 
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Figure  8:  The  dependence  of  scattered  power  (top)  and  the  average  scattering  angle  (bottom)  on  the  number  of 
scattering  objects.  The  width  of  the  simulation  region  is  160  pm,  while  the  length  of  the  simulation  region  is  either 
Lx  =  160  pm  (o)  or  Lx  =  320  pm  (*) 
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Figure  10:  Scattered  light  intensity  from  plane  wave  for  three  sizes  of  cells  that  do  not  contain  a  nucleus :  with 
diameter  Z>0  =  1 1  pm  (dotted  line),  with  diameter  D0  =7  pm  (solid  line)  and  with  diameter  D0  =  5  pm  (dashed  line). 
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Figure  1 1:  Scattered  intensity  from  a  Gaussian  beam  with  f/3  optics  on  a  logarithmic  scale  (a)  and  linear  scale  (b) 
for  three  types  of  cells:  without  a  nucleus  (dotted  line),  with  a  small  nucleus  Dn  —  3  pm  (dashed  line),  and  with  a 
large  nucleus  Dn=6  pm  (solid  line).  In  all  cases,  the  cell  diameter  is  assumed  to  bell  pm  and  the  cell  is  located  in 
center  of  the  focal  spot  of  the  laser  beam.  The  incident  beam  spectrum  is  given  by  the  solid  lines  with  the  *  symbol. 
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Figure  12:  Scattered  intensity  from  a  Gaussian  beam  with  VI 0  on  a  logarithmic  scale  (a)  and  linear  scale  (b)  for 
three  types  of  homogeneous  spherical  cells  of  different  diameters:  with  diameter  Do  —  17  pm  (dotted  line),  diameter 
Do  =  1 1  pm  (dashed  line),  and  diameter  Do  =  8  pm  (solid  line).  In  all  cases,  the  cell  is  located  in  center  of  the  focal 
spot  of  the  laser  beam.  The  incident  beam  spectrum  Is  given  by  the  solid  lines  with  the  *  symbol. 


34 


DRDC  Suffield  TR  2004-187 


DOCUMENT  CONTROL  DATA 

(Security  classification  of  title,  body  oi  abstract  and  indexing  annotation  must  be  entered  when  document  is  classified) 


2.  .  SECURITY  CLASSIFICATION 

(overall  security  classification  of  the  document 
including  special  warning  terms  if  applicable). 


1 .  ORIGINATOR  (the  name  and  address  of  the  organization  preparing  the  document. 
Organizations  for  whom  the  document  was  prepared,  e.g.  Centre  sponsoring  a 
contractor's  report,  or  tasking  agency,  are  entered  in  section  8.) 


Defence  R  &  D  Canada  -  Suffield 

RO.  Box  4000,  Medicine  Hat,  AB,  Canada  T1 A  8K6 


UNCLASSIFIED 


TITLE  (the  complete  document  title  as  indicated  on  the  title  page.  Its  classification  should  be  indicated  by  the  appropriate 
abbreviation  (S,C,R  or  U)  in  parentheses  after  the  title). 

An  Approximate  Method  for  Modelling  Laser  Light  Scattering  from  Biological  Cells 


AUTHORS 

(Last  name,  first  name,  middle  initial.  If  military,  show  rank,  e.g.  Doe,  Maj.  John  E.) 

E.,  Shao,  Y.  and  Yee, 


5.  DATE  OF  PUBLICATION  (month  and  year  of  publication  of  document) 


October  2004 


6a.  NO.  OF  PAGES  (total 

containing  information.  Include 
Annexes,  Appendices,  etc). 


6b.  NO.  OF  REFS  (total  cited  in 
document) 


7.  DESCRIPTIVE  NOTES  (the  category  of  the  document,  e.g.  technical  report,  technical  note  or  memorandum.  If  appropriate,  enter  the  type  of  report, 
e.g.  interim,  progress,  summary,  annual  or  final.  Give  the  inclusive  dates  when  a  specific  reporting  period  is  covered). 

Technical  Report 


8.  SPONSORING  ACTIvrTY  (the  name  of  the  department  project  office  or  laboratory  sponsoring  the  research  and  development.  Include  address). 

Defence  R  &  D  Canada  -  Suffield 

P.O.  Box  4000,  Medicine  Hat,  AB,  Canada  T1 A  8K6 


9a.  PROJECT  OR  GRANT  NO.  (if  appropriate,  the  applicable  research  and 
development  project  or  grant  number  under  which  the  document  was 
written.  Specify  whether  project  or  grant). 

PCN  16QD40 


9b.  CONTRACT  NO.  (if  appropriate,  the  applicable  number  under  which 
the  document  was  written). 


10b.  OTHER  DOCUMENT  NOs.  (Any  other  numbers  which  may  be 
assigned  this  document  either  by  the  originator  or  by  the  sponsor.) 


10a.  ORIGINATOR'S  DOCUMENT  NUMBER  (the  official  document  number 
by  which  the  document  is  identified  by  the  originating  activity.  This 
number  must  be  unique.) 

DRDC  Suffield  TR  2004-187 


1 1 .  DOCUMENT  AVAILABILITY  (any  limitations  on  further  dissemination  of  the  document,  other  than  those  imposed  by  security  classification) 

(X)  Unlimited  distribution 

(  )  Defence  departments  and  defence  contractors;  further  distribution  only  as  approved 
(  )  Defence  departments  and  Canadian  defence  contractors;  further  distribution  only  as  approved 
(  )  Government  departments  and  agencies;  further  distribution  only  as  approved 
(  )  Defence  departments;  further  distribution  only  as  approved 
(  )  Other  (please  specify): 


12.  DOCUMENT  ANNOUNCEMENT  (any  limitation  to  the  bibliographic  announcement  of  this  document.  This  will  normally  correspond  to  the  Document 
Availability  (1 1).  However,  where  further  distribution  beyond  the  audience  specified  in  (1 1)  is  possible,  a  wider  announcement  audience  may  be 
selected). 


UNCLASSIFIED 


1 3  ABSTRACT  (a  brief  and  factual  summary  of  the  document.  It  may  also  appear  elsewhere  in  the  body  of  the  document  Itself.  H  ^  h'9^  d®s^le  flat  ,he 
abstract  of  classified  documents  be  unclassified.  Each  paragraph  of  the  abstract  shall  begin  with  an  indication  of  the  security  c  assfflcaflon  of  the 
information  in  the  paragraph  (unless  the  document  itself  Is  unclassified)  represented  as  (S).  (C),  (R),  or  (U).  It  is  not  necessary  to  include  here  abstracts  In 
both  official  languages  unless  the  text  is  bilingual). 

A  new  technique  for  approximating  the  scattering  of  laser  light  by  biological  cells  is  reported.  This  technique  is 
based  on  a  three-dimensional  scalar  wave  equation  approximation  of  the  full  Maxwell’s  field  equations  for  the 
electromagnetic  field.  This  scalar  wave  equation  for  describing  the  light  scattering  patterns  of  cells  containing 
arbitrary  morphological  structure  (e.g.,  various  organelles)  is  solved  numerically  using  a  spectral  method. 
The  accuracy  of  the  spectral  numerical  method  is  verified  by  comparison  with  solutions  obtained  from  linear 
perturbation  theory  and  Mie  theory.  Comparison  with  Mie  theory  shows  that  the  three-dimensional  scalar 
wave  equation  is  a  good  approximation  to  the  full  Maxwell's  field  equations  for  light  scattering  up  to  moderate 
forward  scattering  angles  (i.e.,  for  scattering  angles  less  than  about  35°).  The  approximate  technique  used 
here  is  capable  of  correctly  predicting  the  scattered  intensity  patterns  from  biological  cells  over  a  dynamic 
range  spanning  six  orders  of  magnitude.  The  new  technique  can  be  applied  to  calculate  the  light  scattering 
either  from  an  individual  biological  cell  and  from  a  sample  containing  an  ensemble  of  such  biological  cells. 
The  scattering  intensity  patterns  predicted  using  the  new  technique  can  potentially  be  applied  to  diagnose 
the  size  and  internal  structure  of  biological  cells,  making  it  a  valuable  interpretative  tool  in  flow  cytometry 
(e.g.,  in  the  detection  of  rare  event  cells  such  as  those  resulting  from  a  biological  warfare  agent  attack,  or  for 
the  rapid  noninvasive  optical  assessment  of  tissue  pathology  in  the  detection  of  cancerous  cells). 


14  KEYWORDS.  DESCRIPTORS  or  IDENTIFIERS  (technically  meaningful  terms  or  short  phrases  that  characterize  a  Oocument  and  couid  bo  h  p 
cataloguing  the  document  They  should  be  selected  so  that  no  security  classification  is  required.  Identifiers,  such  as  equipment  model  designa  lo 
name  ^military  project  code  name,  geographic  location  may  also  be  included.  If  possible  keywords  should  be  selected  from  a  Pu“!^d  ^®““rus.  ®B9. 
Thesaurus  o7Eng!neering  and  Science  Terms  (TEST)  and  that  thesaurus-identified.  If  it  not  possiWe  to  select  indexing  terms  which  are  Unclassified,  the 
classification  of  each  should  be  indicated  as  with  the  title). 

light  scattering 
biological  cells 
spectral  method 
lasers 


UNCLASSIFIED 


Defence  R&D  Canada 


R&D  pour  la  defense  Canada 


Canada's  leader  in  defence  Chef  de  file  au  Canada  enR&D 

and  national  security  R&D  pour  la  defense  et  la  securite  nationale 


www.drdc-rddc.gc.ca 


